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Abstract 

We present a finite element approximation for the one-sided Stefan problem and 
the one-sided Mullins-Sekerka problem, respectively. The problems feature a fully 
anisotropic Gibbs-Thomson law, as well as kinetic undercooling. Our approxima- 
tion, which couples a parametric approximation of the moving boundary with a 
finite element approximation of the bulk quantities, can be shown to satisfy a sta- 
bility bound, and it enjoys very good mesh properties which means that no mesh 
smoothing is necessary in practice. In our numerical computations we concentrate 
on the simulation of snow crystal growth. On choosing realistic physical parameters, 
we are able to produce several distinctive types of snow crystal morphologies. In 
particular, facet breaking in approximately crystalline evolutions can be observed. 
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1 Introduction 

Pattern formation during crystal growth is one of the most fascinating areas in physics 
and materials science. Furthermore, crystallisation is a fundamental phase transition, and 
a good understanding is crucial for many applications. In this paper we will concentrate 
on a mathematical model based on the one-sided Stefan and Mullins-Sekerka problems, 
for which we will introduce a new numerical method of approximation. The numerical 
solutions presented here are tailored for the description of snow crystal growth. However, 
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we note that with minor modifications our approach can be used for other crystal growth 
scenarios; see Barrett et al. (2010b), which in particular have applications in engineering 



as, for example, in the foundry industry. 

The basic mathematical model for crystal growth involves diffusion equations in the 
bulk phases together with complex conditions at the moving boundary, which separates 
the phases. Depending on the application, either heat diffusion or the diffusion of a 
solidifying species has to be considered. If a pure, e.g. metallic, substance solidifies, then 



the basic diffusion equation is the heat equation for the temperature, see Gurtin (1993); 



Barrett et al. (2010b). Whereas for snow crystal growth the diffusion of water molecules 

In the case that a 



Libbrecht (2005). 



in the air is the main diffusion mechanism, see 
binary metallic substance solidifies, then models involving both heat and species diffusion 
simultaneously, and which are coupled through the interface conditions, are considered, 



see e.g. Davis (2001). 



At the moving boundary a conservation law either for the energy or for the matter 
has to hold. In the case of heat diffusion, one has to take into account the release of 
latent heat through the well-known Stefan condition, which relates the velocity of the 
interface to the temperature gradients at the interface; the latter being proportional to 



Gurtin 


(1993 


); 


Davis 


(2001 


); 


Barrett et al. 


(2010b 



growth the continuity equation at the interface relates its velocity to the particle flux at 
the interface, which is given in terms of the gradient of the water molecule density. In 
conclusion, mathematically very similar conditions arise in both models. 

Beside the above discussed continuity equation, another condition has to be specified 
at the interface. In the case that heat diffusion is the main driving force in the bulk, ther- 
modynamical considerations lead to the Gibbs-Thomson law with kinetic undercooling at 
the interface, see Gurtin (1993); Davis (2001); Barrett et al. (2010b). This law relates the 



undercooling (or superheating) at the interface to the curvature and the velocity of the in- 
terface. In the case of snow crystal growth one has to consider a modified Hertz-Knudsen 
formula, which relates the supersaturation of the water molecules at the interface to the 



curvature and velocity of the interface, see e.g. equations (1) and (23) in Libbrecht (2005). 
The physics at the interface depends on the local orientation of the crystal lattice in space, 
and hence the parameters in the interface conditions discussed above are anisotropic. In 
particular, the corresponding surface energy density leads, through variational calculus, 
to an anisotropic version of curvature, which then appears in the moving boundary condi- 
tion, see Giga (2006). In addition, kinetic coefficients in the moving boundary condition 



will also, in general, be anisotropic. 

In the numerical experiments in Section [5j we focus on snow crystal growth, where 
the unknown will be a properly scaled number density of the water molecules. However, 
straightforward modifications, e.g. choosing different anisotropics, allow our approach to 
apply in the context of other crystal growth phenomena. In addition, we note that our 



approach can be used for many other moving boundary problems, see e.g. Barrett et al. 

(polobl. 
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Figure 1: The domain f2 in the case d = 2. 



In earlier work, the present authors introduced a new methodology to approximate 
curvature driven curve and surface evolution, see |Barrett et al\ ( |2007b|a| |2008b[ ). The 
method has the important feature that mesh properties remain good during the evolution. 
In fact, for curves semidiscrete versions of the approach lead to polygonal approximations, 
where the vertices are equally spaced throughout the evolution. This property is impor- 
tant, as most other approaches typically lead to meshes which deteriorate during the 
evolution and often the computation cannot be continued. The approach was first pro- 
posed for isotropic geometric evolution equations, but later the method was generalized 
to anisotropic situations, Barrett et al. ( 2008a[c ) 



geometry was coupled to bulk fields, Barrett et al. (2010b|). 
possible to show stability bounds. In Barrett et al. 



and to situations where an interface 
In most cases it was even 



(2010b) the two-sided Stefan and 



Mullins-Sekerka problems, as a model for dendritic solidification, were numerically stud- 
ied. The physical parameters, such as the heat conductivity, had to be chosen the same 
in both phases; whereas, in this paper we focus on the situation where diffusion can be 
restricted to the liquid or gas phase, respectively. Hence, we need to study a one-sided 
Stefan or Mullins-Sekerka problem. In particular, an anisotropic version of the one-sided 
Mullins-Sekerka problem is relevant for snow crystal growth, 
Section 



see 



Libbrecht (2005) and 



5.1 



below. This, and the fact that the anisotropy in snow crystal growth is so 
strong that nearly facetted shapes occur, makes this application a perfect situation in 
order to test whether our approach is suitable for one-sided models for solidification. 

Before discussing our numerical approach and several phenomena, which we wish to 
simulate, we formulate the anisotropic one-sided Stefan and Mullins-Sekerka problem 
with the Gibbs-Thomson law and kinetic undercooling in detail. Let Q C M d be a given 
domain, where d = 2 or d = 3. We now seek a time dependent interface (r(t)) te r 5?i, 
r(t) CC f2, which for all t e [0,T] separates Q into a domain Q + (t), occupied by the 
liquid/gas, and a domain fi_(i) := Q \ £l+(t), which is occupied by the solid phase. See 
Figure [I] for an illustration. For later use, we assume that {T(t)) te ^ ^ is a sufficiently 
smooth evolving hypersurface parameterized by x(-,t) : T — > IR d , where T C lR d is a given 
reference manifold, i.e. T(t) = x(T,t). Then V := Xt ■ v is the normal velocity of the 



3 



evolving hypersurface T, where v is the unit normal on T(t) pointing into Q + (t). 



We now need to find a time and space dependent function u defined in the liquid/gas 
region such that u(.,t) : Q+(t) —> M and the interface (r(t)) tg [ ^ fulfill the following 
conditions: 



•d u t — /C Au 
du 
du 
pV 

m 

u 

r(o) 



/ 



-AV 



a k 1 — au 



r , 



i? M (.,0) = du 



in fl + (t) , 


(1.1a) 


on r(f) , 


(1.1b) 


on r(t) , 


(1.1c) 


on <9f2 , 


(l.ld) 


in fi+(0) ; 


(Lie) 


possible forcing 


term, while 



T CC Q and Uq : f2+(0) — > R are given initial data. We always assume that the solid 
region 0_ (t) is compactly contained in Q. 

The unknown u is, depending on the application, either a temperature or a suitably 
scaled negative concentration. The orientation dependent function f3 is a kinetic coeffi- 
cient, 7 is the anisotropic surface energy and i? > 0, /C, A, p, a, a > are constants whose 
physical significance is discussed in Section 5.1 and in Barrett et al. (2010b). For snow 
crystal growth, see Section [5TT| — u is a suitably scaled concentration with —ud being the 
scaled supersaturation. 

It now remains to introduce the anisotropic mean curvature x 7 . One obtains x 7 as 
the first variation of an anisotropic interface free energy 



d-l 



7(j?) m 



where 7 : W 1 — > R>o, with j(p) > if p 7^ 0, is the surface free energy density which 
depends on the local orientation of the surface via the normal z7; and 'H d ~ 1 denotes the 
(d — l)-dimensional Hausdorff measure in Mr. The function 7 is assumed to be positively 
homogeneous of degree one, i.e. 



7 (6 p) = 6 7 (p) WpeR d , Wbe 



i{p)-P = l{v) Vpei rf \{0}, 



where 7' is the gradient of 7. The first variation of |T| 7 is given by, see e.g. Giga (2006) 



and Barrett et al. (2008c), 



x 7 := -V S .7V) ; 



where V 5 . is the tangential divergence of T, i.e. we have in particular that 



dt 



mi- 



d 
dt 



j(P) dH 



d-l 



r(t) 



x^VdH 



d-l 



r(t) 



We remark that in the isotropic case we have that 

lip) = lisoip) := \p\ WpeR d 



[1.2) 



;i.3) 
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which implies that j(u) = 1; and so |T| 7 reduces to |T|, the surface area of T. Moreover, in 
the isotropic case the anisotropic mean curvature x 7 reduces to the usual mean curvature, 
i.e. to the sum of the principal curvatures of T. 



In this paper we are interested in anisotropies of the form 

L 



lip) = ^2ht(p)], nip) ■= \P-G t j$ 



;i.4) 



where Ge G 



>dxd 



for i — 1 — y L, are symmetric and positive definite matrices. We note 



that (|1.4|) corresponds to the special choice r = 1 for the class of anisotropies 

' L 



i{p) 



, com 



;i.s) 



which has been considered by the authors in Barrett et al. (2010b). Numerical methods 



based on anisotropies of the form (1.5) have first been considered in Barrett et al. (2008a) 



and Barrett et al. (2008c), and there this choice enabled the authors to introduce un- 
conditionally stable fully discrete finite element approximations for the anisotropic mean 
curvature flow, i.e. (1.1c) with a = 0, and other geometric evolution equations for an 



evolving interface T. Similarly, in Barrett et al. (2010b), the choice of anisotropies (1.5) 



lead to fully discrete approximations of the Stefan problem with very good stability prop- 
erties. We note that the simpler choice r = 1, i.e. when 7 is of the form (1.4), leads to a 



finite element approximation with a linear system to solve at each time level, see (3.6a 



In three space dimensions, the choice (1.4) only gives rise to a relatively small class of 



anisotropies, which is why the authors introduced the more general (1.5) in 



(2008c). For the modelling of snow crystal growth, however, the choice (1.4) is sufficient 



Barrett et al. 



and we will stick to this case in the present paper, but we point out that using the method 



from Barrett et al. (2010b) the approach in this paper can be easily generalized to the 



more general class of anisotropies in (1.5) 



We now give some examples for anisotropies of the form (1.4), which later on will be 



used for the numerical simulations in this paper. For the visualizations we will use the 



Wulff shape, |Wulflj fll901[ ), defined by 

W := {p G R d : p.q < 7(g) VgGM d } 



;i.6) 



Here we recall that the Wulff shape W is known to be the solution of an isoperimetric 



problem, i.e. the boundary of W is the minimizer of 



in the class of all surfaces 



enclosing the same volume, see e.g. Fonseca and Miiller (1991). 



Let L 



[e 2 \p\ 2 + p\ (1 — e 2 )} 2 for e > 0. Then a hexagonal anisotropy in M 2 can 



be modelled with the choice 



lif) = IheM ■= J>(#(0O + f )P) , 



;i.7) 
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Figure 2: Wulff shape in M 2 for (O) with e = 0.01 and 9 = 0. 





Figure 3: Wulff shape in M 3 for (1.8) with e 
with e = 0.01 (right). 



0.01 (left). Wulff shape in R 3 for ( fl~9| 



where R(9) denotes a clockwise rotation through the angle 9, and #0 G [0, |) is a parameter 
that rotates the orientation of the anisotropy in the plane. The Wulff shape of (1.7) for 
e = 0.01 and 9 = is shown in Figure |2j 



In order to define anisotropies of the form (1.4) in 

cos 6» sin6» (A / cos0 

matrices R\(9) := I -sine cose o and R2(9) := I o 1 



o 1 



- sin 9 




we introduce the rotation 



Then 



T (p) = Z £ (i? 2 (f )p) + £ J e (ife(0 o + f ) j3) 



is one such example, where #0 G [0, |) again rotates the anisotropy in the X\ — x 2 plane. 



The anisotropy (1.8) has been used by the authors in their numerical simulations of 



anisotropic geometric evolution equations in Barrett et al. (2008c, 2010c|a ), as well as for 
their dendritic solidification computations in Barrett et al. (2010b). Its Wulff shapes for 
0.01 is shown on the left of Figure |3l 



e = U.Ul is snown on tne lelt ol f igure A small modification of (1.8), which is more 
relevant for the simulation of snow flake growth, is 



lip) = iHeM ■= W(I)P) + 73 i>(J2i(*o + ®& 



;i.9) 



Its Wulff shape for e = 0.01 is shown on the right of Figure [3j We note that the Wulff 



shape of (1.9), in contrast to (1.8), for e — > approaches a prism where every face has the 



same distance from the origin. In other words, for (1.9) the surface energy densities in the 



basal and prismal directions are the same. We remark that if Wq denotes the Wulff shape 



of (1.9) with e = 0, then the authors in Gravner and Griffeath (2009) used the scaled 



Wulff shape \ Wo as the building block in their cellular automata algorithm. In addition, 
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Figure 4: Wulff shape for the approximation of (1.10) with 7 TB = 1 (left) and 7tb 
(right) for e = 10~ 2 . 



= 0.1 



Klett 



we observe that the choice (1.9) agrees well with data reported in e.g. Pruppacher and 



(1997, p. 148), although there the ratio of basal to prismal energy is computed as 



0.92 < 1. 



In addition, we consider an example of (1.4), where L = 2 and G\ = diag(l, l,e 2 
G2 = 7tb diag(e 2 ,5 2 , 1), so that it approximates for small e the anisotropy 



lip) = 7tb N + {p't +P2) 5 



;i.io) 



as considered in e.g. Giga and Rybka (2004). See Figure |4J where we show its Wulff shape 
for 7tb = 1 and 7tb = 0.1 for e = 10 -2 . We note the Wulff shape of (1.10) is given 
by a cylinder with basal radius one and height 27tb- 
diameter is 7tb- 



Hence its ratio of height to basal 



More examples of anisotropics of the form (1.5) can be found in Barrett et al. (2008a|c 



2010c). 



Crystal growth in general, and snow crystal growth in particular, is a highly anisotropic 
mechanism. In snow crystal growth the morphologies that appear depend strongly on 
the environment and, in particular, on the temperature and the supersaturation, which 
influence the values of a and ud, respectively, in (l.la-e). This can be seen in the 



famous Nakaya diagram, see Figure[5j Depending on these parameters, either solid prisms, 
needles, thin plates, hollow columns or dendrites appear in snow crystal growth. The 
anisotropy of the surface energy can be responsible for the hexagonal symmetry, but 
probably also an anisotropic j3 has an influence on the shapes appearing in snow crystal 



growth, see e.g. Libbrecht (2005) and Yokoyama (1993). Depending on the size of the 



crystal, either the kinetic anisotropy or the anisotropy in the surface energy dominates, see 
Yokoyama and Sekerka ( 1992 ) or Kobayashi and Giga (2001 ). It is one of the goals of this 
paper to study the influence of the anisotropics in (3 and 7 on the growth morphologies. It 



was discussed in Libbrecht (2005) that the kinetic coefficient can vary drastically between 
the directions of the two basal hexagonal facets and the directions of the six prismal facets. 
Depending on the environmental conditions either flat crystals or columns crystals appear, 
see Figure [5} 



Gurtin 


(1993 


) and 


Davis 



(2001). The evolution of interfaces driven by anisotropic curvature has been studied by 
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Figure 5: The Nakaya diagram illustrates which snow crystal forms appear at different 



temperatures and supersaturations. This figure is taken from Libbrecht (2005). 



many authors, and we refer to Giga (2006) for an overview. For the full problem (l.la-e), 
to the knowledge of the authors, no existence result seems to be known. Although there 
are results for two-sided variants, see Luckhaus (1990) in the isotropic case and Garcke 



and Schaubeck (2011) in the anisotropic case. We remark that also cases where the Wulff 
shape is crystalline, i.e. has sharp corners and flat parts, have been studied. In this case 



nonlocal curvature quantities have to be considered, and the geometric equation (1.1c) for 



the interface is of singular diffusion type. Giga and Rybka ( |2004 ) showed the existence of 
self-similar solutions for (l.la-e) in a situation, where the Wulff shape is a cylinder. We 



will attempt to compute self-similar solutions in Section [5} 

In snow crystal growth often flat parts appear, and in some cases they become unstable 
and break, see Figure|5 Libbrecht (2005) and Gonda and Yamazaki (1982 ). Only recently, 
have researchers studied facet breaking from a mathematical point of view. In the three 



dimensional case there have been studies by Bellettini et al. (1999) and Giga and Giga 



(1998) for geometrical evolution equations - see also the numerical studies in Barrett 
et al. (2008c). A full crystalline model of solidification facet breaking has, so far, only 



been studied analytically by Giga and Rybka (2006) and numerically by Barrett et al 



(2010b). Clearly from the Nakaya diagram facet breaking is an important issue in snow 



crystal growth, and we will study this aspect numerically in Section [5} 

Numerical approaches for dendritic solidification, that are based on the Stefan prob- 
lem with the Gibbs-Thomson law, are often restricted to two space dimension, see e.g. 



Yokoyama and Sekerka (1992); Schmidt (1998) and Bansch and Schmidt (2000), where in 



the latter article the coupling to a fluid flow is also considered. The first implementations 



in three space dimensions are due to Schmidt (1993, 1996), and the present authors later 



proposed a stable variant of Schmidt's approach which could also handle the anisotropy 



in a more physically rigorous way, see Barrett et al. (2010b). We also would like to re 



fer to the fascinating results on snow crystal growth, which were established by Gravner 
and Griffeath (2009), using a cellular automata model. They were able to compute a 
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large variety of forms, which resemble snow crystals in nature; even though the overall 
approach does not stem from basic physical conservation laws, and it is difficult to relate 
its parameters to physical quantities. 

The outline of this paper is as follows. In Section [2] we introduce a weak formulation 
of the one-sided Stefan problem and the one-sided Mullins-Sekerka problem, that we 
consider in this paper. Based on this weak formulation, we then introduce our numerical 
approximation of these problems in Section [3j In particular, on utilizing techniques from 
Barrett et al. (2010b), we derive a coupled finite element approximation for the interface 
evolution and the diffusion equation in the bulk. Moreover, we show well-posedness 
and stability results for our numerical approximation. Solution methods for the discrete 
equations and implementation issues are discussed in Section |4| In addition, a non- 



dimensionalization of a model for snow crystal growth from Libbrecht (2005 ), which allows 
us to derive physically relevant parameter ranges, is recalled in Section 5.1 Finally, we 



present several numerical experiments, including simulations of snow crystal formations 
in three space dimensions, in Section |5j 



2 Weak formulation 



In this section we state a weak formulation of the problem (l.la-e) and derive a formal 



energy bound. Recall that °& > and /C, A, p, a, a > are physical parameters that are 



discussed in more detail in Section 5.1, below, and in Barrett et al. (2010b) 



We introduce the function spaces 

3 0t+ (t) := {0 G if 1 (£]+(£)) : = on dVl} 
and S D ,+ (t) := {(f) G iJ 1 (fi+(t)) : = u D on dVl} . 

In addition, we define 

V:=H\r,R d ) and W := H l (T, R) , 
where we recall that T is a given reference manifold. A possible weak formulation of 



(l.la-e), which utilizes the novel weak representation of x 7 V introduced in Barrett et al. 



(2008c), is then given as follows. Find time dependent functions u, x and x 7 such that 
u(-,t) G S Dj+ (t), x(-,t) G V, x 7 (-,£) G W and 



+ /C(Vw,V0)+ -(/,$. 



-K. 



Y(t) 



du 
dV 



d-l 



I d-1 



a / x t .v<pm° 

'r(t) 

V0G5 o ,+ (t), (2.1a) 



P 



x t .v% 



d-1 



'r(t) 



T(t) P(H) 

x^p.f) m d - 1 + ( vf x, vf ff). 



[ [ax^- aulxdn^ 1 V X e W , (2.1b) 
Jr(t) 

(2.1c) 



•■'(*) 

o V fjeV 
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hold for almost all times t e (0, T], as well as the initial conditions (Lie). Here (•,•)- 
denotes the L 2 -inner product on Q + (t). 



We note that, for convenience, we have adopted a slight abuse of notation in (2.1a 



c). Here, and throughout, this paper we will identify functions defined on the reference 
manifold T with functions defined on T(t). In particular, we identify v G W with vox^ 1 on 
r(t), where we recall that T(t) = x(T,t), and we denote both functions simply as v. For 
example, x = id is also the identity function on T(t). In addition, we have introduced th e 



shorthand notation (V s -, Vr_) 7 for the inner product defined in |Barrett et al.\ ( |2008c [). 
In particular, on recalling (1.4), we define the symmetric positive definite matrices Gt 
with the associated inner products (•, -)g on ^ d ^ 

Gi := [det G$ \Gi\~ 1 
Then we have that 



and (v,w)q = v . G^w V v , w e 



1 — )■ L 



d-1 



(2.2) 



where 



d-1 



i=i 



with {tj , . . . ,^d-i\ being an orthonormal basis with respect to the Gi inner product for 
the tangent space of T(t); see Barrett et al. (2008c) for further details. 

Assuming, for simplicity, that the Dirichlet data up i s cons tant, we can est ablish the 

following formal a priori bound. Choosing = u — ud in (2.1a), y = -x t . P in (2.1b) and 

- \ -* \ k a ^ 

77 = — in (2.1c) we obtain, on using the identities 



fie) 



d 
dt 



n+(t) 



^ d£ c 



n+(t) 



gVdH 



d-l 



(2.3) 



r(t) 



with C d denoting the Lebesgue measure in M d , see e.g. Deckelnick et al. (2005), and 
d 



dt 



d 
d7 



7 (^)d^- 1 = (Vfx,Vff t ) 



r(t) 



7 ' 



(2.4) 



sec 



Barrett et a/. (2008c), that 



d 
dT 



P - u D\q, + 



a A 



+ 



Xp 



r(t) 



r(t) | 7 - A m d vol(fi+(t)) ) + AC (V u, V u). 
V 2 

- - = -- / 



/3(z7 



/ V| M - MD | 2 d^ d - 1 + (/, M - MD ) + , (2.5) 

JT(t) 



where | ■ |q + denotes the L 2 -norm on £l + (t). In particular, the bound (2.5) for •& > gives 
a formal a priori control on u and T(t) only if V > 0, i.e. when the solid region is not 
shrinking. 
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3 Finite element approximation 



Let = to < t± < . . . < £m-i < *m = T be a partitioning of [0, T] into possibly variable 
time steps r m := t m+ i — t m , m = — > M — 1. We set r := max m= o_>Af-i T m . First we 
introduce standard finite element spaces of piecewise linear functions on Q. 

Let Q be a polyhedral domain. For m > 0, let T" m be a regular partitioning of Q into 
disjoint open simplices, so that Q = \J ™ e q-mo m . Let Jq be the number of elements in 
T m , so that T m = {o™ : j = 1 — » J^ 1 }. Associated with T m is the finite element space 

5 m := { x G C(Q) : x |o»» is linear V o m G T m } C tf 1 ^) . (3.1) 

Let be the number of nodes of T m and let {p t J l }j=i be the coordinates of these nodes. 

K 777 — 

Let {^J 1 } .^ be the standard basis functions for S 1 " 1 . We introduce I m : C(f2) — >■ S 1 ™, 
the interpolation operator, such that (I m f])(p^) = rj(p£) for k = 1 — > K r £. A discrete 
semi-inner product on C(Q) is then defined by 

( Vl , V2 ) h m := (I m [ Vim ],l) , 

with the induced semi-norm given by |ry|n im := [(v^v)m]^ f° r V e C(Q). 

The test and trial spaces for our finite element approximation of the bulk equation 



(2.1a) are then defined by 



S™ ■= { x e S m : X = on dQ} and S£ := { X G 5 m : X = I m u D on dQ} , (3.2) 

where in the definition of we allow for ud G Hz(d£i) fl C(dfi). Without loss of 
generality, let {0™},=!° be the standard basis functions for S™. 



The parametric finite element spaces in order to approximate x and x 7 in (2.1a 



defined as follows. Similarly to Barrett et al. (2008b), we introduce the following discrete 



spaces, based on the seminal paper Dziuk (1991). Let T m C R be a {d — 1) -dimensional 



arc 



polyhedral surface, i.e. a union of non-degenerate [d— l)-simplices with no hanging vertices 
(see Deckelnick et al. (2005, p. 164) for d = 3), approximating the closed surface T(t m ), 
m — — > M. In particular, let T m = {J^aJ 1 , where {cr! f n } J £ 1 is a family of mutually 
disjoint open (d — l)-simplices with vertices {q^ 1 }^- Then for m = — > M — 1, let 

V(T m ) := {x G C(T m ,R d ) : x\ a? is linear V j = 1 -)• Jf n } =: [W{T m )] d C if 1 (r m ,M (i ) , 

where W(T m ) C /f 1 (r m ,M) is the space of scalar continuous piecewise linear functions 
on r m , with {xT}k=i denoting the standard basis of W(T m ). For later purposes, we 
also introduce n m : C(r m ,IR) — > W(T m ), the standard interpolation operator at the 
nodes {q^} k ^, and similarly vf m : C(T m ,R d ) ->■ V(r m ). Throughout this paper, we will 
parameterize the new closed surface r m+1 over T m , with the help of a parameterization 
Xrn+i G ]/(r m ) ; i.e. T m+1 = X m+1 (V m ). Moreover, for m > 0, we will often identify X m 
with id G V_{T m ), the identity function on T m . 
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For scalar and vector functions v , w E L 2 (T m ,R^) we introduce the L 2 inner product 
•) m over the current polyhedral surface T m as follows 



(v,w) m := v.w dU^ 1 . 

Here and throughout this paper, -W denotes an expression with or without the superscript 
*, and similarly for subscripts. If v, w are piecewise continuous, with possible jumps across 
the edges of {crJ l } J £ 1 , we introduce the mass lumped inner product (•, as 

jm , 

w ) h m ■= 3 E K\ X> • «o((ffi~)> ( 3 - 3 ) 
i=i ^=1 

where {<^™}fc =1 are the vertices of aj 1 , and where we define t>((g™)~) := lim 

Here laj 1 ] = mzjv \(QjZ ~ Qjl) A • ■ • A (q^ — q^)\ is the measure of aj 1 , where A is the 
standard wedge product on R d . Moreover, we set | ■ 1^^) := (•, -)m • 

Given T" 1 , we let Q™ denote the exterior of T m and let fi™ denote the interior of T m , 
so that T m = dVL™ = fi™ fl fFT. In addition, we define the piecewise constant unit normal 
v™ to T m by 

(fin _ gm\ A . . . A ftfm _ gm\ 



-7 



K^-^)A---A(g--g-) 



where we have assumed that the vertices {^}^ =1 of a™ are ordered such that v" 1 : T" 1 — >■ 
M d induces an orientation on T m , and such that zT" 1 points into 



Before we can introduce our approximation to (2.1a-c), we have to introduce the notion 
of a vertex normal on T m . We will combine this definition with a natural assumption that 
is needed in order to show existence and uniqueness, where applicable, for the introduced 
finite element approximation. 

(A) We assume for m = -> M - 1 that \af\ > for all j = 1 -»■ JfT\ and that T m C H. 
For fc = 1 -> if™, let := {a™ : $™ G a™} and set 



A 



-= U^es"^. and w fc .= — ^ I j I v i 



Then we further assume that u™ ^ 0, = 1 — > K™, and that dimspanju;™}^ = d, 
m = -> M - 1. 



Given the above definitions, we also introduce the piecewise linear vertex normal 
function 



k=l 
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and note that 

(v,wu m )^= (v,wLJ m )^ ^ veV(T m ), w eW{T m ). (3.4) 



Following Barrett and Elliott (1982), we consider the following unfitted finite element 



approximation of (2.1a-c). First we need to introduce the appropriate discrete trial 



and test function spaces. To this end, let be an approximation to Q™ and set 

Qrn,h ._ ^ ^ {}™' h . We stress that need n 
T m . Then we define the finite element spaces 



Qrn,h ._ ^ ^ {}™' h . We stress that fi™ ,?t need not necessarily be a union of elements from 



rltti. 



S£ := { X £S m : X (pT) = if supp0f C 



nm cm p. om 



(3.5) 



Our finite element approximation is then given as follows. Let T , an approximation 
to T(0), and, if ■& > 0, U° G S° D be given. For m = -»■ M - 1, find £/ m+1 G S] 
Xm+i E y(Y m ) and G VL(r m ) such that for all <p G x e VL(r m ), 77 G Z( r 



■m 



)C(VU m+ \V<p) 



X ( 7T r 



m,+ 1 



P 



-1* 



m+1 



O]" 1 - ^,X^ m ) -«« +1 ,x) m + a(f/ m+1 , X ) 



ill , 



and set T m+1 = x m+1 (r m ). Here we define 



(3.6a) 
(3.6b) 
(3.6c) 



(Vx,V^) m , + := / V X .V^d£ d = V^^ / V X .V^d£ d V X ,^S 

O m ' h 3 In™ 

JU + 3=1 Jo j 



(3.7a) 



and, in a similar fashion, 



ixMt, + :=E^|SP / /m [x^]d/: d V Xl ^S™. (3.7b) 

3=1 j 



For later use we note that it follows immediately from (3.5) and (3.7a) that 



(V^,V^) m ,+ >0 V<^GS m + \{0} 



(3. 



In addition, we set f m+1 (-) '■— f(-,t m +i), where we assume for convenience that / is 
defined on Q. In addition, for $ > 0, U° G Sp is given by U° = I°[u Q ], w here u G C(Q) 
is an appropriately defined extension to Q of the given initial data from (Lie). 
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Moreover, (V s •, V s -) 7)Tn in (3.6c) is the discrete inner product defined by 



(3.9) 



Note that (3.9) is a natural discrete analogue of (2.2), see Barrett et al. (2008c) for details. 



This choice of discretization will lead to unconditionally stable approximations in certain 



situations; see Theorem 3.1, below. 



Remark. 3.1. We note that for $ > the approximation (3.6a-c) is only meaningful 



when the discrete solid region does not shrink. To see this, assume that the discrete solid 
region shrinks at some time step, so that S™, \ Sq 1 ^ 1 ^ for some m > 1. Assume for 
simplicity that T m = T m_1 , so that S m = S m -\ Now let <j>f e \ S^ 1 , which means 
that the node p™ is an active node in S™ + , but was inactive in S™+ , i.e. U m (pJ L ) = 
since U m E S^^ 1 ■ Here the value U m (pJ 1 ') = is arbitrary, and has no physical meaning. 
Cr uciall y, however, this value will play a role on the discrete level, since choosing <p = (fff 1 
in (3.6a) ; and noting that ((ft™ , </>™)m + > ®> means that JJ m+1 will depend on U m (pj l ). 



In practice this technical restriction is not very relevant, since in physically meaningful 
simulations the solid region typically never shrinks. Here we also recall that the formal 



energy bound (2.5), for d > 0, is also only meaningful, when the solid region is not 
shrinking. 

Theorem. 3.1. Let the assumption (A) hold. Then there exists a unique solution 
(U m+1 ,X m+1 ,n™ +1 ) e x V(T m ) x W(T m ) to ( |3.6af c). Let u D e R and define 

cmljjm v-ra\ . ^ \TTm |2 Ot \ i_ m i 

t [U ,A ).--\U - u D \ nm+ + — |1 | 7 , 



(3.10) 



where I. 



f2,m, 



' /m,+J 



. Then the solution to (3.6a-c) satisfies 



cm/Trm+l \?m+]\ i \ „, / \?m+l \?m ,-tm\h i u lrrm+1 jjm\2 

t [U ,A ) + Au D {A -A ,U ) m + -\U —U \^ m _ 



+ T m K(yU m+ \VU m+1 ) m , + + r, 
<E m (U m ,X m )+r m {f 



m + 2 

Xp ' 
a 

m+l jjm-\-l 



if 71 )} 



1 X 

2 



m+l 



x r * 



T„ 



m,h 



(3.11) 



Proof. As the system (3.6a-c) is linear, existence follows from uniqueness. In order to 



establish the latter, we consider the system: Find (U,X,Kj) 6 S™ + x V_(T m ) x W(T r 
such that 



'm,+ 
TV. 



# (V, + r m K (V U, V <p) m>+ -A {n r 

ii 



X.uj" 



V y? G S m + , (3.12a) 



a{K l7 xr m + a(U, X }m = V X G W(F 



(« 7 co m , ff) h m + (Vf X, Vf ff) Jtm = V ff E m r 



(3.12b) 
(3.12c) 
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on noting (3.4), that 



Choosing ip_^U in (|3.12a[), x = ^ m [X . u m ] in ( |3.12b| ) and 77 = ^ X in ( |3.12c| ) yields, 

V [/3(z7")]~iX.u; m +^(V^X,V s G X) 7 , m 

(3.13) 



(U, U)*+ + r m K (V U, V + 



a A 



0. 



It immediately follows from (3.13) and (3.8) that U = € 5™ 



. In addition, on recalling 



that a, A > 0, it holds that X = X c e R d . Together with (|3. 13 ), for U = 0, and the 
assumption (A) this immediately yields that X = 0, while (3.12b) with x = k 7 implies 
that k 7 = 0, compare Theorem 3.1 in Barrett et a/. (2008c). Hence there exists a unique 



solution (U m+ \X 



m+l \^m+l ^.m+l 



e x v(r m ) x H/(r m ). 



It remains to establish the bound (3.11). Let denote the characteristic function of 



a set A. Choosing <p = U m+l - u D I m AL™, h in (3.6a), X = ~ a vr m [(X m+1 - X m ) . Q m ] in 



( 3.6b) and 77 = ^ (X m+1 - X m ) in ( |3.6c| yields that 
(f/^+i _ [/™ £/™+i _ + r m /C (V £/ m+1 , V t/ m+1 ; 

+ — (V? X m+1 , V? (X m+1 - X m )) 7 , m + r m ^ 
a a 

*m+l 



m,- 



m,h 



\ „. I vm+1 vm ,- t .m\h 1 _ / x 

-A«o (A ^ -X ,a; ) TO + r m (/ 



em+l jjm+l \h 



and hence (3.11) follows immediately, where we have used the result that 



(vfx m+ \ vf (x m+1 - x m )) 7 , m > |r 



m+l 



I pm I 

7 \*- \y 1 



see e.g. Barrett et al. (2008a) andlBarrett et al. (2008c) for the proofs for d = 2 and d = 3, 



respectively. □ 

The above theorem allows us to prove unconditional stability for our scheme under 
certain conditions. 



Theorem. 3.2. Let the assumptions of Theorem \3.1\ hold with up = 0. In addition, 
assume that 1 
it holds that 



assume that either *& = or that U m G S™ and Vl™' h C f2+ 1,h for m = 1 — > M — 1. Then 



£ m {U m+1 ,X m+1 ) + J2 T k£ 



k=0 



(V U k+1 ,VU k+1 ] 



k,+ + 



Xp 



1 x k+l - x k 



Tk 



k,h 



< 



£°(U ,X ) + J2r k (f k+1 ,U k +%, 



(3.14) 



fc=0 



for m — — >■ M — 1 . 
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Proof. The result immediately follows from (3.11) on noting that, if i? > 0, it follows 
from our assumptions that £ m (U m ,X m ) < £ m ~ 1 (U m ,X m ) for m — 1 — > M — 1, since 
then 



i m [(u m y] dc d < 



i m [(u m y] dc 



rm—l 



m\2-\ 



dC c 



(3.15) 



□ 



Remark. 3.2. Theorem 3.2 establishes the unconditional stability of our scheme (3.6a-c) 
under certain conditions. Of course, if ud ^ ; analogous weaker stability results based 
on (3.11) can be derived. We note that the condition U m G is trivially satisfied if 



S 



rn—T 
D 



c S 



D , e.g. when mesh refinement routines without coarsening are employed. 



m—l,h 
+ 



The 

on the other hand, is ensured whenever the discrete solid region 



condition Q™' h c Q 

is not shrinking. This is in line with the corresponding continuous energy law (2.5). 
Note also that the condition U m G S™, would be too strong, as in physically meaningful 
computations the solid region grows, and so the condition would enforce that U m = at 
vertices which are now in the solid region, but were degrees of freedom in S^l . In the 



simpler case that $ = 0, the stability bound (3.11) is independent of U m and so here the 



stability bound (3.14) holds for arbitrary choices of bulk meshes T r ' 



REMARK. 3.3. With the techniques introduced in this paper, it is a simple matter to 



extend the finite element approximation introduced in Barrett et al. (2010b) for the two- 
sided Stefan problem with constant heat conductivity fC = K s = /Q to the case fC s — fCi ^ 0; 



where we have adopted the notation from Barrett et al. (2010b, (2.1a-e)). Here the 



subscripts s and I refer to the solid and liquid phase, respectively. 

Our finite element approximation for this problem is then given as follows. Let T be 
given. For m = ->■ M - I, find U m+1 G S r g, X m+1 G V(T m ) and k™ +1 G W(T m ) such 
that for all ip e S m , x G W(T m ), ff G V(T m ) 



Jjm+l _ jjr. 



T,, 



+ 



i£{l,s} 



[^(vc/- +i ,v^) mji -(/r +i ^)y 



m+l 



x r 



TV? 



7T 



x m+l - x r 



. Id 



a(^ +1 ,x) h m + a(U m+1 , X ), 



„m+l ,-tm 



m ,ff) h m + (V?X m+1 ,V?ff) 1 , m = 0. 



(3.16a) 

(3.16b) 
(3.16c) 



an 



dsetT m+l =x m+1 (r m ). Here (Vx,V^) m ,j and (x, v) h m i; for i G {s,l} and for x^ e 



S m , are defined analogously to (3.7a,b), where Q, 



■\in,h 



{l™> h and n™> h :-- 



represent 



approximations to the "liquid" and "solid" phases in this two-sided Stefan problem. 
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4 Solution of the discrete system 



Introducing the obvious abuse of notation, the linear system (3.6a-c) can be formulated 

\ ( U m+1 \ 



as: Find (U m+1 , 5X m+1 ) such that 



f^Mn + An 

m 





where (U m+1 , 





a M r 



[N T 



Tm L T 

At 



J 



K 



m+1 
m+1 



\SX m+1 J 



( i M n U m + o m \ 




(4.1) 



-1 5 Xm+l) G 



O x 



X 



o^if™ denote coefficients of these 



finite element functions with respect to the standard bases of S m , W(T m ) and V^_(T m ), 
respectively. The definitions of the matrices in (4.1) directly follow from (3.6a-c), but we 
state them here for completeness. Let i, j 



1 ->■ Kg and k, I 



1 -> K^. Then 



[Mn] ii :=tf(C,Ct,+ > 
[M r , Q ]« := (C,xP>m, 

[Nr)ki ■= (xT,xT^ m ) h m , 



[A 



n\ij 



/C(V0f,V. 



((c^[(xre-;-).^])r 



1 < i < ^H!d 



m\ -jrr 



rifci •- 



[ivf ]„ : = m^T'xTiXu^Tm = m^VxT.XkYm^l 



(Vf (xrel),Vf ( X re,)) 7 , 



»J=1 



l„,m „,m\h 



(4.2) 



where {e*j}f =1 denotes the standard basis in IR d and where we have used the convention 
that the subscripts in the matrix notations refer to the test and trial domains, respectively. 
A single subscript is used where the two domains are the same. We note that the special 
definition of An, together with g m in (4.1), accounts for the Dirichlet boundary conditions 
of U m+1 G Here g m is defined by 



9i 




1<*<K% D , 
K% D <i< KR 



(4.3) 



Clearly, the matrix Aq will in general be singular. In particular, it will have zero diagonal 
entries for every vertex pj 1 e fi™' . Hence we enforce U' m+1 G + by setting 



[A 



n\ij 



[A 



5. 



1,3 



[Anh^O, 
[A n }ij = ; 



(4.4) 



i.e. we replace zero diagonal entries by 1. 



The assembly of the matrices in (4.2), apart from Aq, is described in 



Barrett et al. 



(2010b, §4). The assembly of Aq, and in particular the possible d efini tions of the region 



n 



will be discussed in Section 



4.1 



below. The linear system (4.1) can be efficiently 



solved with iterative solvers applied to a Schur complement formulation, see Barrett et al. 



(2010b) for details. For completeness we state that for the application of preconditioners 
and for the solution of subproblems we make use of the packages LDL and AMD, see 



Davis (2005); Amestoy et al. (2004). 
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4.1 Definition of the discrete liquid/gas region 



We now discuss possible choices of in (3.5) and (3.7ab). To this end, we partition 



the elements of the bulk mesh T m into liquid/gas, solid and interfacial elements as follows. 
Let 



' + 


:= {o m 


e r m 


-m 


en;}, 


'j-m 


:= {o rn 


e T m 


-m 


c n m } , 


<-i—m 

1 pm 


:= {o m 


e r m 


-m 


n r"V 0} 



(4.5) 



Then T m = T+ U T m U 7f™ is a disjoint partition. 



Clearly, using ap n = is not very practical, since the intersection of f2™ with 
elements of the bulk mesh T m can be complicated. Moreover, computing the domain Q™ 
is unlikely to be rewarded with lower overall approximation errors, since the trial and 



test functions in (3.7ab) are only piecewise linears. Instead, we consider the following 
approach, which defines Vt™' h with the help of a piecewise linear approximation to 

as 

Q™> h : = {p E Q : {I m X^){p) > 0} . (4.6) 



Next we discuss an algorithm that computes for the strategy (|4.6|). We assign 



each element of T m to one of the three sets T) 



+ ' 



T m or J$k as follows: 



Algorithm 1: Mark all bulk mesh elements as liquid/gas, solid or cut. 



find all elements of 7^™ . 
and 77" := 0. 



1. Traversing over T r 

2. Set T:=T m \ 7~ r " m 

3. Move all elements that touch dfl from T to ■ 

4. For as long as this is possible, move neighbours of elements in 

5. Set T_ m := r. 



from T to T\ 



+ • 



In addition, for later use, we need to decide for each bulk mesh vertex p" 2 
whether it belongs to ff™ or to This can be done as follows. 



3 ^ 



7TL 777, 

Algorithm 2: Assign all bulk mesh vertices to fl_ or Q + . 

1. All vertices of elements in T™ belong to Q_. 

2. All vertices of elements in belong to Q + . 

3. For any remaining vertices {pj 1 }, choose a pj 1 with a neighbouring vertex q that 
is known to belong to Vt_ or Vt + . If T m cuts the segment [pT 1 , q] C M an even 
number of times, assign ftp to the same region as q; otherwise to the opposite 
region. Repeat this, until all vertices have been assigned. 

Remark. 4.1. The global Algorithm [T] is only needed at the very first time step. For sub- 
sequent time steps, the existing marking of bulk mesh elements can be updated depending 
on the movement ofT m . In particular, only elements in 7^7 1 \7p^ need to be considered. 



18 



On assuming that T m has not travelled over a whole bulk mesh element, these elements 
can be marked with the help of neighbour information. This is far more efficient than 
employing the global Algorithm [T] at every time step. 

In addition, in practice for a refined bulk mesh in the neighbourhood of T m , all the 
remaining vertices in Step 3 of Algorithm [2] have immediately a neighbouring vertex q 
that is known to belong to Q_ or Q + . 



An alternative approach to (4.6) would not define Q™ ,h explicitly, but rather the effect 
of Q™' h on the inner products defined in (3.7a,b). Her e it is natural to define Q™' h in such 
a way, that {J m eT mO m C Vl^ h . Then the integral in (3.7a) can be rewritten as 



v(o m ) / Vx-Vy?d£ d 



(4.7) 



where v(o m ) G [0,1] denotes the fraction of the element o m that is consider ed to belong 
to the liquid/gas region Q™' h ; and similarly for the inner product defined in (3.7b). Note 
that (4.7) only implicitly defines (candidates of) the region 



In practice, several choices of v(o m ) G [0,1] can be considered. The approach (4.6) 
corresponds to 

v (o m ) = 1 V o m G T r Z , (4.8a) 



while the choice 



v (o m ) = V o m G 



(4.8b) 



was used in Bansch and Schmidt (2000) for a two-sided Stefan problem with nonvanishing 



heat conductivity coefficients. We note that for the one-sided situation considered in this 



paper, the strategy (4.8b) does not make sense, as it dramatically affects the accuracy of 



the approximation JJ m+1 on T m . An alternative approach is the choice 



v{o m ) 



k 
d+1 



I m dC d 



V O m G Tr m 



(4.8c) 



where k denotes the number of vertices of o m that lie within Q + . A simpler approach is 
to set 



v(o m ) 



V O m G 7r 



111 



(4.8d) 



We note that for the practical implementation, the strategies (4.8a b,d) only need the 
marking from Algorithm [TJ The additional Algorithm [2] is only required for the strategy 
(4.8c). In practice, the three strategies (4.8a, c,d) all show very similar numerical results. 



Hence, in general we will employ the simplest strategy (4.8a) 
Remark. 4.2. Of course, setting 



v ( m ) = \o m n sr_ 



V o m G 7?£ 



(4.9) 
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corresponds to Vt™ ,h = f2™. This will in general be too costly to do in practice. However, 
we mention one possible strategy here. For an arbitrary open bounded set V C M. d with 
Lipschitz boundary it holds that 



vol(V) 



1 dC c 



(id — zq) . V v dH 



d-l 



(4.10) 



dv 



z G 



where id is the identity function on I 

Vy denotes the outer normal to V. Applying (4.10) for V 

m p| pm 



d is an arbitrarily fixed point, and where 
H Q™, on noting that 



Vy 



on o 



and Vy = v m , the outer normal of o r ' 



on do m n tt m 



yields a 



way of using (4.9) in practice. Of course, in th is case V is a polytope, with dV being a 
union of flat facets. Thus the integral in (4.10) simplifies on noting that id . Vy is now 
constant on each facet, and vanishes on each facet that contains z . Moreover, o m fl T m 
can be computed as in Barrett et al. (2010b, §4-5)- 



where for our purposes it is enough to compute |-F„nf2+| for \i 



It remains to calculate do m fl Vt™, 
1 — > d + 1, where F„ are 



the edges/faces of o r ' 



i.e. 



do r ' 



USF,. Ford 



of F^ n which is straightforward. For d = 
disjoint union of possibly non-convex polygons. 



2 this reduces to finding the lengths 
3 the set F^ fl Q™ in general can be the 
The oriented boundary of these polygons 



can be found by suitably arranging the line segments making up dFn fl Y m , as well as the 



line segments making up dF^ fl Q T : 
Gauss ' area formula. 



Then the area 



Ffj, n ST 



can be easily computed with 



5 Numerical results 



We implemented our finite element approximation (3.6a-c) within the framework of the 
finite element toolbox ALBERTA, see Schmidt and Siebert (2005). We use the bulk 
mesh and parametric mesh refinement strategies introduced in Barrett et al. (2010b, §5). 



Here the bulk mesh adaptation algorithm, was inspired by a similar strategy proposed in 
Barrett et al. (2004) and Banas and Niirnberg (2008) for d = 2 and d = 3, respectively, 
results in a fine mesh of uniform mesh size hf around T m and a coarse mesh of uniform 
mesh size h c further away from it. Here hf = and h c = ^ are given by two integer 

numbers Nf > N c , where we assume from now on that Q = (-H, H) d . For the one-sided 



problem s considered in this paper, we slightly amend the strategy from Barrett et al. 
(2010b, §5), in that we allow an even coarser grid inside Q m,h . Of course, the definitions 



(3.5) mean that this has no effect on the numerical results. Moreover, the parametric 



mesh refinement uses bisections in order to avoid elements getting too large over time. 
We stress that apart from this simple mesh refinement, no other changes were performed 
on the parametric mesh in any of our simulations. In particular, no mesh smoothing 
(redistribution) was required. 



r, m 



Throughout this section we use (almost) uniform time steps; in that, r m 
— > M — 2, and Tm-i = T — t m -i < t. Unless otherwise stated we set Q = 
with H = 4. Similarly, unless otherwise stated, we always employ the strategy (|4.8al) for 



-H,HY 



the computation of VL™' h . The initial interface T(0) is always a circle/sphere of radius 
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Rq G (0,if) around the origin. For the Stefan problem, i.e. if $ > 0, we set 



«oW = 

unless a true solution k is given. 
For later purposes, we define 





X _ e Ro-\A 
1 _ e Ro-H 
u D 



\z\ < R , 
ud Ro < \z\ < H , 
\z\>H; 



X(t) :- 
and similarly for U . 



t—tm-X j£ m I t m —t Y-m—1 



Tm-X 



x r 



(5.1) 



t e [t 



m— 1; "m\ 



m > 1 



5.1 Non-dimensionalization of a model for snow crystal growth 



An aim of this paper is to be able to perform computations for the growth of snow crystals 
with realistic parameters and on physically relevant length and times scales. Upon non- 



dimensionalizing the continuum model for snow crystal growth from Libbrecht (2005), it 



turns out that (l.la-c) with 

= 0, K=l, A = l, p = 1.42x10" 



a 



10" 



1, / = o 



(5.2) 



is a physically realistic model. Here the typical length scale is 100 /xm, typical time scales 
vary from 100 s to 1300 s, — u denotes a scaled concentration of water vapour in the gas 
phase and — ud is a scaled supersaturation. We refer to Barrett et al. (2012) for more 
details on the physical interpretation of these parameters. 



5.2 Convergence experiments 



We begin with a compa rison of the approximation error vol(f2 + (0)) — vol(Q+ h ) for the 

Q \ Bi(0) and, for the case d = 2, 
2 7+% and iV>. = 4*. An example of 



four different strategies (4.8a-d). Here we set f2+(0) 
use the spatial discretization parameters Nf = - 
how the discrete interface T cuts the bulk mesh T° is shown in Figure |6j The numerical 
results are shown in Table [TJ where we observe that the strategies (4.8c d) produce far 
smaller errors than (4.8a b). However, as we will see in the subsequent convergence 
experiments, this does not seem to have an influence on the overall approximation error 
for the underlying solutions u and V . For completeness, we repeat the same experiments 
for d = 3, where now N f = 2 6+i , iV c = 4* and K° = K(i), with (K(0), K(l), K(2), K(3)) = 
(770,3074, 12290,49154), for i = 3. The results are shown in Table g 
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Figure 6: Parts of the triangulation T° and the interface T when Nf 

— ll 2 

> 2) ■ 



N c = 4. From left to right [-2,2] 2 , [-1,0] 2 and [-1 



2 8 and 



h 



f 



6.2500e-02 
3.1250e-02 
1.5625e-02 
7.8125e-03 
3.9062e-03 



hf 



5.4874e-02 
2.7439e-02 
1.3720e-02 
6.8600e-03 
3.4300e-03 



(4.8a) 



-2.3534e-01 
-1.1425e-01 
-5.5655e-02 
-2.7579e-02 
-1.4273e-02 



jl8c) 



-1.1384e-02 
-4.8739e-03 
-1.9442e-03 
-1.2118e-03 
-7.0317e-04 



(4.8d) 



Table 1:0= (-4,4) 2 . Approximation error vol(O+(0)) - vo\(Q^ h ) for (4.8a -d). 



-8.7802e-03 
-4.8739e-03 
-1.4559e-03 
-7.2351e-04 
-8.7610e-04 

0,h\ 



(4.8b) 



2.1778e-01 
1.0450e-01 
5.2743e-02 
2.6132e-02 
1.2521e-02 



5.2.1 One-sided Stefan problem 



Next we investigate the approximative properties of our algorithm (3.6a-c) for the fol- 



lowing exact solution to the one-sided Stefan problem ( l.laj -e), in the case of an isotropic 
surface energy, i.e. 7 = | ■ |. Here we adapt the following expanding circle/sphere solution 



for the two-phase Stefan problem in Barrett et al. (2010b, (6.5)), where the radius of the 
circle/sphere is given by r(t); and so fl+(t) = fl\ B r ^(0). Assume that 



and let 



r{t) 



tf = /C = X = p 



! (o) + ^ 



a 



w(t) 



r(t) 



v{s) = - 



1 

C4 



s -I? 2 

e 4 z 



dz 



h 



f 



hf 



(4.8a) 



(4.8c) 



(4.8d) 



(4.8b) 



1.2500e-01 
6.2500e-02 
3.1250e-02 
1.5625e-02 



2.0854e-01 
1.0472e-01 
5.2416e-02 
2.6215e-02 



-8.9192e-01 
-4.5246e-01 
-2.2370e-01 
-1.1247e-01 



-6.7696e-02 
-1.7403e-02 
-3.5485e-03 
-8.0954e-04 



-1.2902e-03 
-3.2433e-03 
4.1878e-04 
-2.8311e-04 

0,h\ 



8.8933e-01 
4.4598e-01 
2.2454e-01 
1.1190e-01 



Table 2: Vt = (-4,4) 3 . Approximation error vol(f2+(0)) - vol(O^) for (4.8a -d). 
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hf 

J 


1 


\\U -I h u\\ La °, + 


| \X — x\\l°° 




1 


6.2500e-02 


5.0640e-02 


2.4595e-01 


9.2545e-02 


677 


128 


3.1250e-02 


2.7093e-02 


7.2888e-02 


2.1049e-02 


1329 


256 


1.5625e-02 


1.3740e-02 


2.0818e-02 


3.5439e-03 


2753 


512 


7.8125e-03 


6.8637e-03 


5.2596e-03 


6.2892e-04 


8853 


1024 


3.9062e-03 


3.4307e-03 


1.2318e-03 


2.1081e-04 


71305 


2048 



Table 3: Q = (-4,4) 2 and T = 1. Convergence test for (5.3) with (4.8a). 





(4.8c 


) 


(4.8d 


) 


hf 


\\U - I h u\\ L oo !+ 


=f 

\\X — 2 1 Loo 


||[/-/ ft «|| LO c i+ 


=f 

| \X — x\\l°° 


6.2500e-02 
3.1250e-02 
1.5625e-02 
7.8125e-03 
3.9062e-03 


2.5105e-01 
7.7931e-02 
2.3909e-02 
6.1309e-03 
1.7662e-03 


9.9965e-02 
2.4780e-02 
4.8305e-03 
1.5131e-03 
7.6027e-04 


2.5204e-01 
7.8798e-02 
2.4363e-02 
6.2823e-03 
1.8835e-03 


1.0185e-01 
2.5502e-02 
5.1527e-03 
1.7784e-03 
9.4426e-04 



Table 4: Q = (-4,4) 2 and T = 1. Convergence test for (5.3) with (4.8c) and (4.8d) 



Then it is easy to see that on letting 



dt 



w{t) 



2r 3 (t) 



the solution u to (l.la-e), with ud in (l.ld) replaced by u\g D n, is given by the restriction 
to n+(t) of 

'w{t) ?6fl.(i), 

w{t) > - (5 ' 3) 



u(z, t) 



r(t) 



For d 

where we use r(0) = Ro 
4* and r = 4 s "* x 10~ 3 . 



2, we perform the following convergence experiment for the solution (5.3), 

2 7 +\ 



= 0.5. For i = 
The errors \\U 



2K° 



-)> 4, we set N f 

and \\X — x\\l 



I h u\\ L o 



val [0,T] with T = 1, 
J h d| X oo . := max m= i_> M 



so that r(T) « 1.12, are displayed in Table |3 
\\U m -I m - 1 u(t m )\\ 00 ^ + , where ||[/™-/-- ! 



on the inter- 
Here \\U — 



max. 



feAf m-i \U m (p) - u(t m ,p)\ and AT 



m— 1 



| oo,m— 1,+ ) 



Moreover \\X 



max fe= i^ A -s 



mm 



max m= i^ M \\X 

A < -x(y,t m ) 



■m 



x(-, t 



^(^m) || oo,m— 1,+ 
m—X 



j = i -> k™- 1 } no; nn+(t 



where ||X(t T 



and /ip := max - 



diam(cr 



M\ 
5 i 



. Note that 



2 Ky due to the growth of the interface. In addition, we use the convergence 
experiment in order to compare the different strategies (4.8c) and (4.8d). See T able |4| 
where we present the same computations as in Table |3l but now for (4.8c) and (4.8d). 



For the new results we omit the additional mesh statistics, as they are very similar to the 
results for (4.8a) shown in Table [3j 
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J 


± 


\\U -I h u\\ L ™ 


— * 

\\X — X\\l°° 


K 

i L 


1 


6.2500e-02 


5.0474e-02 


2.4940e-01 


9.7039e-02 


645 


128 


3.1250e-02 


2.7082e-02 


7.3208e-02 


2.2291e-02 


1353 


256 


1.5625e-02 


1.3739e-02 


2.0678e-02 


3.9277e-03 


2753 


512 


7.8125e-03 


6.8641e-03 


4.9403e-03 


7.2470e-04 


9017 


1024 


3.9062e-03 


3.4309e-03 


1.2377e-03 


2.8003e-04 


74589 


2048 



Table 5: Q = (— 4, 4) 2 and T — 1. Convergence test for the two phase Stefan problem. 



hf 




\\U-I h u\\ L ^,+ 


\\X — x\\l°° 


K M 




1.2500e-01 


1.1309e-01 


1.9195e-01 


5.1473e-02 


1655 


770 


6.2500e-02 


5.9856e-02 


8.7871e-02 


2.0037e-02 


5353 


3074 


3.1250e-02 


3.0712e-02 


2.8850e-02 


5.2297e-03 


26221 


12290 


1.5625e-02 


1.5464e-02 


8.3717e-03 


1.0781e-03 


356903 


49154 



Table 6: Q, = (-4,4) 3 and T = 0.1. Convergence test for (5.3) with (4.8a). 



in 



We also compare the numbers in Tables [3] and [4] with the corresponding errors for the 
approximation from Barrett et al. (2010b) for the two-phase Stefan problem, see (2.1a-e) 
Barrett et al. (2010b), with the same choice of parameters. Note that u(-,t) : f2 — > R 

The corresponding errors, where 

jjm jrn—l 



as defined in ( J5.3[ ) then is the desired true solution 



\U 



'U\\L° 



max m= i_i.M 



u 



loo , can be seen in Table 5 



Similarly to Table |3| we perform a convergence test for the solution (5.3) to the 
one-sided Stefan problem, now for d = 3, leaving all the remaining parameters fixed as 
before. To this end, for % = -»■ 3, we set N f = 2 6+i , iV c = 4\ K$ = K(i), where 
(K(0),K(1),K(2),K(3)) = (770,3074,12290,49154), and r = A 3 ~ { x HT 3 . The errors 
\\U — I u\\l°°,+ and \\X — x\\l°° on the interval [0,T] with T = 0.1, so that r(T) ~ 0.59, 
are displayed in Table [6j In addition, we use the convergence experiment in order to 
compare the different strategies (4.8c) and (4.8d). See Table[7| where we present the same 
computations as in Table [6l but now for (4.8c) and (4.8d). We also compare the numbers 



in Tables [6] and [7] with the corresponding errors for the approximation from Barrett 
et al. (2010b) for the two-phase Stefan problem with the same choice of parameters. The 



corresponding errors can be seen in Table pi 





(4.8c 


) 


(4.8d 


) 


hf 


\\U-I h u\\ L « + 


=f 

\\X - X\\ L 00 


\\U-I h u\\ L ~ + 


=f 

\\X — X || ^oo 


1.2500e-01 
6.2500e-02 
3.1250e-02 
1.5625e-02 


1 .9608e-01 
8.7603e-02 
2.8999e-02 
9.3255e-03 


5.1569e-02 
2.0783e-02 
5.9048e-03 
1.4821e-03 


1.9699e-01 
9.0314e-02 
3.0329e-02 
9.9485e-03 


5.1769e-02 
2.1104e-02 
6.1388e-03 
1.6143e-03 



Table 7: Q = (-4,4) 3 and T = 0.1. Convergence test for (5.3) with (4.8c) and (4.8d) 
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Table 8: 



hf 


hf 


\\U - I h u\\ L °° 


— * 

\\X — x\\l°° 


K M 




1.2500e-01 


1.1297e-01 


1.9491e-01 


5.2057e-02 




1781 


770 


6.2500e-02 


5.9798e-02 


8.3255e-02 


2.0582e-02 




5353 


3074 


3.1250e-02 


3.0700e-02 


2.7380e-02 


5.4506e-03 


26221 


12290 


1.5625e-02 


1.5462e-02 


8.1295e-03 


1.1521e-03 


356909 


49154 


: Q = (— 4, 4) 3 and T = 0.1. Convergence 


test for the two phase Stefan p 


h f 


hf 


\\U-I h u\\ LB o r 


i- 


|| X — x Hi, 00 


K M 




6.2500e-02 


8.5583e-02 


5.9751e-02 




1.1650e-02 


1005 


128 


3.1250e-02 


4.2909e-02 


3.7601e-02 




1.6311e-02 


1981 


256 


1.5625e-02 


2.1304e-02 


9.0157e-03 




4.0322e-03 


4069 


512 


7.8125e-03 


1.0632e-02 


1.5531e-03 




6.7227e-04 


11149 


1024 


3.9062e-03 


5.3145e-03 


4.7394e-04 




2.0761e-04 


70733 


2048 



Table 9: Q = (-4,4) 2 and T = 1. Convergence test for (5.4) with (4.8a). 



5.2.2 One-sided Mullins— Sekerka problem 



We start with a comparison of our algorithm (3.6a-c) for the following exact solution to 
the one-sided Mullins-Sekerka problem (l.la-e) with $ = 0, in the case of an isotropic 



surface energy. Here we use the following expanding circle/sphere solution, where the 
radius of the circle/sphere is given by r(t). Assume that 



= 0, /C = A = p 



a 



f = 



and let 



r (t) = (r 2 (0) + 2ty 



Then it is easy to see that the solution u to (l.la-e), with ud in (l.ld) replaced by u \a D n, 
is given by the restriction to Q+(t) of 



u(z, t) 



d 

>(t) 

-ln^- 

rW _ 1 _ 



2 
T{t) 

3 
r(t) 



d 
d 



2 
3 



(5.4) 



For d — 2, we perform the following convergence experiment for the solution (5.4), 
where we use r(0) — Rq — 1. For i = — > 4, we set Nf = K° = 2 7+ \ N c = A 1 and 
r = 4 2 ^ 1 x 10~ 3 . The errors \\U — I h u\\l°°,+ and JJX — af||£oo on the interval [0, T] with 
T = 1, so that r(T) w 1.73, are displayed in Table M In addition, we use the convergence 



experiment in order to compare the different strategies (4.8c) and (4.8d). See Table 10 



where we present the same computations as in Table [9j but now for (4.8c) and (4.8d). 



We also compare the numbers in Tables [9] and 10 with the corresponding errors for 
the approximation from Barrett et al. (2010b) for the two-phase Mullins-Sekerka problem 
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(4.8c 


) 


(4.8d 


) 


hf 


\\U-I h u\\ L oo >+ 


=f 

\\X — x\\l°° 


\\U- I h u\\ L oo >+ 


=f 

\\X — x||t<x> 


6.2500e-02 
3.1250e-02 
1.5625e-02 
7.8125e-03 
3.9062e-03 


7.0732e-02 
4.1221e-02 
1.1504e-02 
3.2383e-03 
1.2919e-03 


6.1554e-03 
1.3540e-02 
2.6877e-03 
3.8846e-05 
1.4815e-04 


7.5587e-02 
4.3588e-02 
1.2409e-02 
3.4367e-03 
1.3623e-03 


5.0174e-03 
1.2923e-02 
2.3901e-03 
1.7735e-04 
2.1749e-04 



Table 10: Q, = (-4,4) 2 and T = 1. Convergence test for (5.4) with (4.8c) and (4.8d) 



hf 




\\U -I h u\\ L °° 


X — x \\loo 


K M 




6.2500e-02 


8.5582e-02 


5.5854e-02 


1.1640e-02 


1005 


128 


3.1250e-02 


4.2910e-02 


3.3181e-02 


1.6328e-02 


1981 


256 


1.5625e-02 


2.1305e-02 


8.6904e-03 


4.0428e-03 


4073 


512 


7.8125e-03 


1.0632e-02 


1.5719e-03 


6.8315e-04 


11493 


1024 


3.9062e-03 


5.3145e-03 


4.7787e-04 


2.1309e-04 


79197 


2048 



Table 11: tt 
problem. 



[— 4, 4) 2 and T — 1. Convergence test for the two phase Mullins-Sekerka 



with the same choice of parameters, when the function u(-,t) : O — > R from (5.4) is the 



desired true solution. The corresponding errors can be seen in Table 11 



Similarly to Table |9| we perform a convergence experiment for the true solution (5.4) 
to the one-sided Mullins-Sekerka problem, now for d = 3, leaving all the remaining 
parameters fixed as before. To this end, for % = — > 3, we set Nf = 2 5+l , N c = 4 J , = 
K(i), where (K(0),K(1),K(2),K(3)) = (770,3074,12290,49154), and r = 4 3 ~* x KT 3 . 
The errors \\U — I h u\\l°°.+ and \\X — x\\l°° on the interval [0,T] with T = 0.1, so that 



r(T) w 1.1 are displayed in Table 12 In addition, we use the convergence experiment in 



order to compare the different strategies (4.8c) and (4.8d). See Table 13, where we present 



the same computations as in Table 12, but now for (4.8c) and (4.8d). We also compare 



the numbers in Tables 12 and 13 with the corresponding errors for the approximation 



from Barrett et al. (2010b) for the two-phase Mullins-Sekerka problem with the same 



choice of parameters. The corresponding errors can be seen in Table [L4 



hf 




\\U - I h u\\ La °, + 


|| X — x|| XjOo 


K M 




2.5000e-01 


2.2637e-01 


1.8264e-01 


1.3621e-02 


1437 


770 


1.2500e-01 


1.1441e-01 


8.2741e-02 


2.6208e-03 


4769 


3074 


6.2500e-02 


5.7328e-02 


3.2617e-02 


8.0637e-04 


22659 


12290 


3.1250e-02 


2.8688e-02 


5.8383e-03 


2.4496e-04 


339431 


49154 



Table 12: Q = (-4,4) 3 and T = 0.1. Convergence test for (5.4) with (4.8a) 
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(4.8c 


) 


(4.8d 


) 


hf 


\\U-I h u\\ L c B>+ 


=t 

\\X — x\\l°° 


\\U - I h u\\ L ^ t+ 


=f 

| \X — x\\l°° 


2.5000e-01 
1.2500e-01 
6.2500e-02 
3.1250e-02 


1.7194e-01 
7.1850e-02 
2.9357e-02 
9.6310e-03 


1.5249e-02 
2.3731e-03 
5.3446e-04 
2.7820e-04 


1.7596e-01 
7.8187e-02 
3.2027e-02 
1.0533e-02 


1.4567e-02 
2.8742e-03 
8.1515e-04 
4.1290e-04 



Table 13: Q = (-4,4) 3 and T = 0.1. Convergence test for (5.4) with (4.8c) and (4.8d) 



hf 




\\U-I h u\\ L o, 


\\X - x\\ L oo 


K M 




2.5000e-01 


2.2681e-01 


1.8285e-01 


1.2023e-02 


1563 


770 


1.2500e-01 


1.1458e-01 


6.7414e-02 


1.3748e-03 


4847 


3074 


6.2500e-02 


5.7385e-02 


2.2704e-02 


7.5695e-04 


22773 


12290 


3.1250e-02 


2.8688e-02 


6.4026e-03 


2.5641e-04 


340087 


49154 



Table 14: ft 
problem. 



-4, 4) 3 and T = 0.1. Convergence test for the two phase Mullins-Sekerka 



What all of the numerical results in Tables [3-14 reveal is that the three strategies 



(4.8a c,d) all behave very similar in practice, with the simple strategy (4.8a) surprisingly 
showing the smallest errors in general. Combined with the fact that implementing this 
strategy requires the fewest computational steps, means that from now on we will always 



use (4.8a) in our experiments. 



5.3 Crystal growth simulations for d = 2 



Throughout this subsection we use the parameters in (5.2) and 7 = jhex defined by (1.7) 
with e = 0.01 and 6> = ^. We use this rotation of the anisotropy jhex, so that the 
dominant growth directions are not exactly aligned with the underlying finite element 
meshes T m . For the kinetic coefficient we usually set — 7. Moreover, the radius of the 
initial crystal seed T(0) is always chosen to be Ro = 0.05. 

We begin with a value of up = —0.004. The results are shown in Figure [7j We also 
show the same experiment for (3 = 1, see Figure [8] We observe that for this experiment, 
the kinetic coefficient (3 appears to have hardly any influence on the growth of the crys- 
tal. Moreover, we can observe that the initially circular crystal seed almost immediately 
assumes a shape that is favoured by the anisotropy 7, i.e. a shape that is close to the 
Wulff shape. This shape then expands at first in a self-similar fashion, before dendritic 
arms start to grow at the vertices of the shape. In order to underline the different effects 
of 7 and (3, we compare the results in Figure [8] with an experiment where we reverse the 



choices of 7 and (3; i.e. we choose an isotropic surface energy density 7 = 7j SO as in (1.3) 



while the kinetic coefficient is defined by (3 = 'jhex, recall (1.7). The numerical results for 
this experiment can be seen in Figure 19} 
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Figure 7: (ft = (-4, 4) 2 , u D = -0.004, 7 = = lhex ) X(t) for t = 0, 0.5, . . . , 5 (left), for 
t — 0, 5, . . . , 50 (middle), and for t — 0, 50, . . . , 500 (right). Parameters are iV/ = 256, 
N c = 4, K° = 16 and r = 0.1. 




Figure 8: (ft = (-4,4) 2 , u D = -0.004, 7 = 7?iea; , /3 = 1) X(t) for t = 0, 0.5, ... ,5 (left), 
for t = 0, 5, . . . , 50 (middle), and for t = 0, 50, . . . , 500 (right). Parameters are Nf = 256, 
N c = 4, K° = 16 and r = 0.1. 



Before we look at experiments with larger values of \uo\, we present the results for a 
run with ud = —0.004, but now run on the larger domain ft = (—8, 8) 2 and until the later 
time T = 2500. See Figure 10 for the results, where the different effects of 7 and /3 are once 
again visible. In fact, the results for the isotropic surface energy 7 = ^ iso seem to indicate 
that the orientation of the underlying finite element mesh has a larger influence on the 
directions, in which the unstable interface grows, than the kinetic coefficient /3 = 'jhex 
itself. To confirm this interpretation, we present a further comparison. This time, we 
choose all coefficients as isotropic, so that 7 = 7j SO and (3 = 1. The corresponding result 
is shown on the right of Figure 10 Once again it appears that the role that /3 plays here 
is insignificant. We observe that in the case that 7 is isotropic a tip splitting instability 




Figure 9: (ft = (-4, 4) 2 , u D = -0.004, 7 = 7iso , (3 = lhex ) X(t) for t = 0, 0.5, . . . , 5 (left), 
for t = 0, 5, . . . , 50 (middle), and for t = 0, 50, . . . , 500 (right). Parameters are Nf = 256, 
N c = 4, K° = 16 and r = 0.1. 
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Figure 10: (ft = (-8, 8) 2 , u D = -0.004, 7 = lhex , f3 = 1 (left), 7 = 7iso , /3 = lhex (middle), 
7 = 7iso , /5 = 1, (right)) X(t) for t = 0, 100, . . . , 2500. Parameters are N f = 512, N c = 8, 
K° = 16 and r = 0.1. 




Figure 11: (ft = (-4, 4) 2 , w D = -0.01, 7 = /3 = 7^) X(t) for t = 0, 5, . . . , 50 (left), and 
for t = 0, 50, . . . , 200 (right). Parameters are N f = 512, iV c = K$ = 16 and r = 5 x 10" 3 . 



occurs. 



In the next experiment, we set ud = —0.01 for 7 = = 'jhex- The results are shown in 
Figure 11 and we observe that a larger supersaturation enhances the unstable behaviour. 



In the next experiment, we set «p = —0.04. The results are shown in Figure 12 The 
distribution of U at time t = 40 can be seen in Figure 13 Here we note that, according 
to the definitions (3.5), in these plots U is set to zero inside the solid phase. 



As a comparison, we repeat the same experiment as in Figure 13 now for (i) the one- 




Figure 12: (ft = (-4, 4) 2 , u D = -0.04, 7 = p = j hex ) X(t) for t = 0, 0.5, . . . , 5 (left), and 
for t = 0, 5, . . . , 40 (right). Parameters are Nf = 1024, N c = K® = 64 and r = 2.5 x 10~ 3 . 
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Figure 13: (One-sided Mullins-Sekerka problem) X(t) for for t = 0, 5, ... ,40 (left). X(t) 
and U(t) for t = 40 on [-4,4] 2 (middle) and on [-2,-1] x [-1.1,-0.1] (right). 




o 




Figure 14: (One-sided Stefan problem) X(t) for for t — 0, 5, . . 
for t = 40 on [-4,4] 2 (middle) and on [-2,-1] x [-1.1,-0.1] 



.,40 (left). X(t) and U(t) 
(right). 



sided Stefan problem, (ii) the two-sided Mullins-Sekerka problem, and (hi) the two-sided 
Stefan problem with •§ — 1 for the Stefan problems. Note that for (ii) and (hi) we employ 
the finite element approximation from Barrett et al. (2010b), while for (i) we use (3.6a -c) 
with •§ — 1. The corresponding plots are shown in Figures [TlfTlBj We observe that the 



difference between the one-sided and the two-sided problems is not very pronounced, but 
one notices that the sidearms in the two-sided problems grow slower due to the fact that 
diffusion into the crystal is possible. 

In the hnal experiments for d = 2, we return to the one-sided Mullins-Sekerka problem 
and set ud = —0.08 and ud = —0.2. The results are shown in Figures 17 and 
respectively. 



5.4 Crystal growth simulations for d 



Throughout this subsection, unless otherwise stated, we use the parameters in (5.2) and 
jhex defined by (1.9) with e = 0.01 and 9q 



7 



12' 



Once again, we use this rotation of 
the anisotropy 7, so that the dominant growth directions are not exactly aligned with the 
X\- and X2-directions of the underlying finite element meshes T m . Moreover, the radius 
of the initial crystal seed T(0) is always chosen to be Rq = 0.05. For later use, we define 
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Figure 15: (Two-sided Mullins-Sekerka problem) X(t) for for t — 0, 5, ... ,40 (left). X(t) 
and U(t) for t = 40 on [-4,4] 2 (middle) and on [-2,-1] x [-1.1,-0.1] (right). 




Figure 16: (Two-sided Stefan problem) X(t) for for t = 0, 5, ... ,40 (left). X(t) and U(t) 
for t = 40 on [-4,4] 2 (middle) and on [-2,-1] x [-1.1,-0.1] (right). 




Figure 17: (fi = (-4,4) 2 , w D = -0.08, 7 = /3) X(t) for t = 0, 0.2, ... ,2 (left), and for 
t = 0, 2, . . . , 20 (right). Parameters are N f = 1024, N c = K^ = 64 and r = 10" 3 . 



Figure 18: (Q = (-4,4) 2 , u D = -0.2, 7 = /3) X(t) for t = 0, 0.04, ... ,0.4 (left), and for 
t = 0, 0.4, . . . , 6.4 (right). Parameters are N f = 2048, N c = K® = 128 and r = 2.5 x 10~ 4 . 
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Figure 19: (fi = (-4,4) 3 , u D = -0.004, (3 = 1) X(T) for T = 50. Parameters are 
N f = 128, iV c = 16, = 98 and r = HT 1 . 




Figure 20: (fi = (-4,4) 3 , u D = -0.004, /3 = /3 flat , 3 ) A(T) for T = 50. Parameters are 
N f = 128, N c = 16, K° = 98 and r = HT 1 . 



the kinetic coefficients 

/5flat(p) = /?flat,^(p) := 

and 



10 



with leN 



(5.5a) 



MP) = PtaXiM ■= [l^ 2t {p{+pi)+Pp with left. (5.5b) 

We note that in practice, similarly to the two-dimensional results in Figures [7] and |8j there 
was hardly any difference between the numerical results for a kinetic coefficient (3 that is 
isotropic in the X\ — a^-plane, such as /3fl at and /3 ta u, and one that is anisotropically aligned 
to the surface energy density, such as e.g. f3 = /3fl at 7. Hence in all our experiments we 
always choose coefficients that are isotropic in the x\ — X2-plane, e.g. (5.5a b). 



In the first experiment, we set «o = 
coefficients (3 = 1 and (3 = /3fl a t,3; see Figures fT9| and [20 



0.004 and compare the results for the two 
We observe that the kinetic 



coefficient seems to be responsible for the fact whether solid prisms or thin plates grow, 
see also the Nakaya diagram in Figure [5] and Libbrecht (2005). More details of the 
evolution for the evolution in Figure [20] are given in Figure [2TJ A continuation of the 
evolution shown in Figure 21, now on the larger domain Q = (— 8,8) 3 , can be seen in 



Figure 22 , where the onset of dendritic growth can be observed. 



An experiment for u e> 
prism grows. 



-0.002 and (3 = (3 ta \\ 1 can be seen in Figure 23, where a solid 



An experiment for up = —0.008 and (3 = /5taii.2 can be seen in Figure 24 In this 



case the basal facets break leading to hollow columns, see Figure [5] and Giga and Rybka 
(|2006l. 
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Figure 21: (fi = (-4,4) 3 , u D = -0.004, /3 = /3 flat)3 ) for t = 1, 2, 10, 20, 30, 40, 50. 
Parameters are N f = 128, N c = 16, K% = 98 and r = 10 -1 . 




Figure 22: (Q = (-8, 8) 3 , u D = -0.004, = /3 flat>3 ) X(t) for t = 50, 100, 150, 200. 
Parameters are N f = 256, N c = 32, K$ = 98 and r ='l0 -1 . 



..■■if 




Figure 23: (Q = (-4,4) 3 , u D = -0.002, (3 = Aaii.i) X{t) fort = 1, 2, 5, 10, 20, 30, 40, 50; 
and X(50) within Q. Parameters are N f = 128, N c = 16, K$ = 98 and r = 10 _1 . 




Figure 24: (fi = (-4,4) 3 , w D = -0.008, /3 = /3 t aii, 2 ) X(t) for t = 1, 2, 5, 10, 20, 30, 40, 50; 
and X(50) within fi. Parameters are N f = 128, N c = 16, K° = 98 and r = 10 -1 . 
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Figure 25: (SI = (-4,4) 3 , u D = -0.02, (3 = /3 flat , 3 ) X(t) for t = 0.05, 0.1, 0.2, 0.3; and 
X(0.3) within Q. Parameters are N f = 512, N c = 32, K% = 1538 and r = 5 x 10~ 4 . 




An experiment for uo = —0.02 and (3 = /3fl a t,3 can be seen in Figure 25 In this case 
the prism facets break leading to capped columns which also can be observed in nature, 
see 



Libbrecht (2005). 



An experiment for Ud — —0.02 and /3 = 
7G0 = 2Z e (J?a(f)ji) 



Aiat,3> but for the anisotropy 7 defined by 

3 

^UMOo + ^P). (5.6) 



with e = 0.01 and 9 = ^ can be seen in Figure |26| This leads to a geometrically more 



complicated breaking of the prismal facets. These can also be observed in nature and 
they are called hollow plates, see Libbrecht (2005). 

We also performed simulations varying (3 in time. This is realistic as a growing snow 
crystal falls to the earth through changing weather conditions, which influences the gov- 
erning parameters, e.g. via the temperature. In the first such example, we choose 



m 

In a second example we choose 



Aall,3(p) 



t e [0,30), 
t e [30, 00) . 



(5.7a) 



(5.7b) 



Atat, 3 (p) *€[0,20), 

/W($ *e[20,oo). 

Results for these choices of (3 and for up = —0.004 can be seen in Figure 27 The shapes 
in Figure [27] can also be observed in nature and they are called scrolls on plates. 
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Figure 27: (Q = (-4,4) 3 , u D = -0.004, /3 as in flSJab) ) X{T) for T = 50. Parameters 



are N f = 128, N r 



16, ^ 



98 and r = 10" 



I I 



mn 




Figure 28: {Q = (-4,4) 3 , u D = -0.004, /3 = Aaii,i) for* = 1, 2, 5, 10, 20, 30, 40, 50; 
and X(50) within fi. Parameters are N f = 128, A^ c = 16, K% = 98 and r = 10 _1 . 



The remaining numerical experiments are for the cylindrical anisotropy (1.10) with 
recall Figure [4j The first case is for 7tb = 1, «o = —0.004 and (3 = /3 t , 



10" 



and the results, which show facet breaking both in the basal and prismal directions, can 



be seen in Figure 28 Some plots of the concentration are shown in Figures 29 and 30 



where Berg's effect, see e.g. Giga and Rybka (2003), can clearly be seen; i.e. U increases 



towards the centre of the basal face before facet breaking occurs. 



For the anisotropy (1.10) it is of interest to find for what value of 7tb the evolution 
of (l.la-e) with 



# = 0, K 



A 



P 



a 



= 7, / = 



(5.* 




Figure 29: (fi = (-4,4) 3 , u D = -0.004, (3 = Aaii,i) X{t) n {z : z x = 0} and U{t) | 21=0 for 
t = 15, 25, 50. Parameters are N f = 128, A^ c = 16, K$ = 98 and r = 10 _1 . 
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Figure 30: (Q = (-4,4) 3 , u D = -0.004, (3 = Aaii,i) X{t) D {z : Z X = 0} and C/(t) | Z1=0 for 
t = 10, 15. Parameters are N f = 128, N c = 16, = 98 and r = 10 -1 . 




Figure 31: (Q = (-8,8) 3 , 7 TB = 0.925) X(t) for t = 0, 0.1, 0.2; and X(0.2) within Q. 
Parameters are X/ = 512, iV c = 32, K™ = 1538 and r = 10~ 4 . 



is self-similar. For example, in Giga and Rybka (2004) it was shown that there exists a 



value 7tb > for which this is the case. Numerically this can be checked by starting this 
flow with a scaled Wulff shape (or a shape close to that), and then to observe whether 
the height to basal diameter ratio of the evolving approximate cylinder converges to 7tb- 

In practice we choose T(0) to be a cylinder with basal radius R = 0.1 and a height/basal 
diameter ratio of 7tb- I n order to obtain the desired sign for V, i.e. for an expanding 
evolution, we set Ud = — 21 in (l.ld). For the domain Q we choose Q = (— 8,8) 3 . 



In practice we appear to obtain a value for self-similarity for some 7tb £ [0.92,0.93], 
although the precise value seems to depend on the resolution of the bulk mesh. In Fig- 
31 we plot some results for an experiment with 7tb = 0.925, while in Figure 32 



ure 



we show the evolution of the ratio of interest for two experiments with 7tb = 0.92 and 
7tb = 0.925, respectively. These results seem to indicate that there exists a value 7tb 
close to 7tb — 0.92 for which the evolution of (l.la-e) with (5.8) and (1.10) is self-similar. 
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Figure 32: Plots of the height/basal diameter ratio for the two runs with 7 TB = 0.92 (top) 
and 7tb = 0.925 (bottom). The dashed lines show the value of 7tb- 
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Conclusions 



We have presented a fully practical finite element approximation for one-sided Mullins- 
Sekerka and Stefan problems with anisotropic Gibbs-Thomson law and kinetic undercool- 
ing. In particular, the method allows the approximation of a continuum model for snow 
crystal growth, which is based on rigorous thermodynamical principles and balance laws. 
To our knowledge, the numerical results presented in this paper are the first simulations 
of snow crystal growth that are based on such a rigorous, physically motivated model. 

In our numerical simulations of snow crystal growth in three space dimensions, we 
were able to produce a significant number of different types of snow crystals. In partic- 
ular, recall Figure [5j we obtained results that resemble solid plates, solid prisms, hollow 
columns, dendrites, capped columns and scrolls on plates. Also facet breaking in the 
moving boundary problems computed have been observed in cases with nearly crystalline 



anisotropic energies, see also Giga and Rybka (2006) for theoretical predictions of facet 
breaking. We therefore believe that the results presented here may help to understand 
the different factors that play a role in the shaping of snow crystals in the real world. 

Producing more complicated dendritic shapes in three space dimensions, with compli- 



cated substructures such as steps and ridges, as in e.g. Libbrecht (2005, Fig. 1), or as in 



the beautiful simulations in Gravner and Griffeath (2009), which were obtained with a cel- 



lular automata algorithm, would need a much higher computational cost when computed 
with the help of a discretized moving boundary problem for a diffusion equation. The 
main reason is that the highly detailed and irregularly structured surface of snow flakes, 
see e.g. Figure 1(c) in Libbrecht (2005), would need to be accurately captured with a 
triangulated surface T m , say. On this surface, a second order partial differential equation 
then needs to be solved, which is coupled to a PDE in the bulk. The necessary resolutions 
for both meshes, as well as the involved computational effort to solve the linear systems 
arising from (3.6a-c), mean that on currently available computer hardware those kind of 



computations cannot be performed. 

Nevertheless, it is our belief that the numerical methods presented here, combined 
with suitable randomizations and fluctuations of physical parameters, together with so- 
phisticated computing equipment should be able to produce all the possible variations of 
realistic snow crystals. In addition, we believe that the computations presented in this 
paper are the most accurate and complex which have been computed so far with the help 
of a Stefan or Mullins-Sekerka problem with hexagonal symmetry. 

Acknowledgement: We are thankful to Prof. Libbrecht for allowing us to use Figure [5} 
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